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In the recent work of S. Sharma et al., (arxiv.org: cond-matt/0912.1118), a single-electron spec- 
trum associated with the natural orbitals was defined as the derivative of the total energy with 
respect to the occupation numbers at half filling for the orbital of interest. This idea reproduces the 
bands of various periodic systems using the appropriate functional quite accurately. In the present 
work we apply this approximation to the calculation of the ionization potentials and electron affini- 
ties of molecular systems using various functionals within the reduced density-matrix functional 
theory. We demonstrate that this approximation is very successful in general and in particular for 
certain functionals it performs better than the direct determination of the ionization potentials and 
electron affinities through the calculation of positive and negative ions respectively. The reason for 
this is identified to be the inaccuracy that arises from different handling of the open- and closed-shell 
systems. 

PACS numbers: 31.15.vo 71.15.-m 



I. INTRODUCTION 

It is generally accepted today that the Fermi surfaces of 
metallic systems obtained with density functional theory 
(DFT), even at the level of local density approximation 
(LDA), arc in good agreement with experiments. Un- 
fortunately, this is not the case with the band gaps of 
insulators and semiconductors which are highly underes- 
timated by most of the exchange-correlation (xc) func- 
tionals within DFT. Even with the exact xc functional 
of DFT, the Kohn-Sham (KS) gap is not expected to 
reproduce the experimental gap[l|. This deviation from 
experiment is most dramatic for Mott-insulators, most of 
which arc predicted by their KS spectrum to be metallic 
while they are experimentally known to be insulating in 
nature. 

In this regard reduced density matrix functional the- 
ory (RDMFT) has shown great promise in improving on 
DFT results for a wide class of systems in that it not 
only improves the KS-band gaps for insulators in gen- 
eral, but also predicts the correct insulating nature for 
Mott insulators 0. Within RDMFT the total energy of 
a system of interacting electrons is expressed in terms of 
the one-body reduced density matrix (1-RDM), 7(1-, r'). 
This energy functional is then minimized with respect 
to 7 under the 7V-representability conditions Q which 
restrict the minimization to the domain of 1-RDMs that 
correspond to ensembles of A^-electron wave functions. A 
major advantage of RDMFT comes from the fact that the 
exact kinetic energy is easily expressed as a functional of 
the 1-RDM of the ground state. In addition, due to the 
departure from the idempotent single-determinant solu- 
tion, static electronic correlations are well described [1]. 
The total ground-state energy as a functional of 7 reads 



(atomic units are used throughout) 
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£[7] =-7iJ J-™ V r7(r, r') dV + J p(r) U cxt (r) d 3 r 
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where p(r) = 7(1% r). U ext is a given external potential, 
and E xc we call the xc energy functional. In practice, the 
xc functional is an unknown functional of the 1-RDM and 
needs to be approximated. A milestone in the develop- 
ment of approximate functionals of the 1-RDM is the 
Miillcr functional [1, which has the following form: 



E xc [ 7 }=E xc [{<j> j },{n j }} = 

If /-|7 1/2 (r,r')| 2 
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where 1/2 is an exponent in the operator sense. Diagonal- 
ization of 7 produces a set of natural orbitals (the eigen- 
vectors of 7), <pj, and occupation numbers (the eigen- 
values of 7), nj. The Miiller functional is known to 
over correlate [6l-il0l]. however, there exist several other 
approximations most of which are modifications of this 
functional and are known to improve results for finite 
systems gEMl- 

Several of these RDMFT functionals reproduce the 
discontinuity of the chemical potential at integer num- 
ber of electrons which is a measure of the fundamental 
gap of the system 1(|_3Q|, UK- More precisely, it was 
demonstrated fl6l. T3C 7131 1 that the complete removal of 
the self-interaction (SI) terms leads to discontinuities in 
the chemical potential that are in good agreement with 
the fundamental gap for finite systems. Unfortunately, 
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this removal has no effect in the total energy for infinitely 
extended natural orbitals in periodic systems since their 
contribution vanishes in the limit of the size of system 
going to infinity. To overcome this problem, Sharma et 
al. @ introduced the power functional [H, [32[ that repro- 
duces discontinuities without requiring the removal of SI 
terms. This functional has the form 

E^[ 1 \=E^[{<j> i },{n j }] = 

where a is an exponent in the operator sense. The power 
functional was applied in the calculation of the funda- 
mental gap of various systems 0, [27j including transition 
metal oxides 0). An optimal value of a between 0.6 and 
0.7 was found to reproduce gaps of all systems in close 
agreement with experiments. These gaps were obtained 
from the discontinuity of the chemical potential, fi(N), at 
an integer total number of electrons N. A problem of this 
method of predicting the gap is that the shape of (Jt(N) 
differs substantially from a step function leading to large 
error-bars in the prediction of the gap. A second prob- 
lem is that one needs to calculate the total energy (and 
//) for several values of N making the calculation time 
consuming. Finally, a third problem is that this method 
does not allow for the calculation of quantities other than 
the gap, for direct comparison with experiments, like for 
example the density of states for extended systems and 
ionization potentials (IPs) and electron affinities (EAs) 
for finite systems. 

An advantage of DFT is that the KS eigenvalues can 
be used as an approximate single-electron spectrum of 
the system. Thus, quantities like the IP and EA can be 
easily estimated using the KS spectrum. A fundamental 
difference between RDMFT and DFT is the lack of a KS 
system within RDMFT and the lack of eigenvalue equa- 
tion makes it difficult to obtain (even approximately) the 
IPs and the EAs. One way to calculate IPs in RDMFT, 
is to use extended Koopman's theorem (EKT), as was 
proposed by Pernal and Cioslowski(33j. They demon- 
strated that the Lagrangian matrix in RDMFT is identi- 
cal with the generalized Fock matrix entering EKT. Thus, 
IPs can be calculated by diagonalization of this matrix. 
They used this idea in the calculation of IPs for small 
molecular systems using the so called Buijse Baerends 
Corrected (BBC) @ and Goedecker Umrigar (GU)[H1 
functionals and showed that the error in the obtained IPs 
is of the order of 4-6%. The same idea was employed in 
combination with yet another xc functional, namely the 
PNOFlfTHUi} functional, for the calculation of the first 
IPs (FIPs) as well as higher IPs (HIPs) of molecular sys- 
tems yielding results of similar quality [34j. However, the 
application of this method is restricted to finite systems, 
since, for solids, it would require the diagonalization of a 
large matrix in wave- vector space. 

Sharma et al. in Ref. (35j introduced an alternative 
technique to obtain spectral information. We refer to this 



Funct. 


A DEF (%) 


A D if(%) 


A DE r(%) 


Muller[5j 


13.24 


12.82 


10.03 


GU[11] 


5.90 


10.19 


9.19 


Power [2] 


12.51 


8.68 


6.08 


AC3[22J 


5.98 


9.33 


6.33 


PNOFl[13,JJ 


6.27 


11.05 


7.21 


BBC3[8] 


6.24 


10.19 


8.79 


ML [19] 


6.15 


9.22 


4.17 



TABLE I: Average absolute errors Adef, Adif, and Ader, 
in the calculation of IP, FIPs and HIPs, for a set of atoms 
and molecules with calculations performed using various xc 
functionals in conjunction with DEF, DIF, DER methods re- 
spectively. 

technique as the "derivative" (DER) method as it entails 
for each natural orbital, k, the associated energy, efc, be 
obtained as the derivative of the total energy with respect 
to the occupation number, rife, at = 1/2 and with the 
rest of the occupation numbers set equal to their ground- 
state optimal values. This technique has been applied for 
the calculation of densities of states of transition-metal 
oxides (NiO, FeO, CoO and MnO) and the results were 
found to be in excellent agreement with experiments (35j 
and other state-of-the-art many-body techniques like dy- 
namical mean-field theory and the GW method. 

In the present work this technique is applied to finite 
systems. We discuss the validity of the approximations 
necessary for the accuracy of the method. We present re- 
sults for the FIP as well as HIP, for atoms and molecules 
as well as the EAs of atoms, molecules, and radicals 
adopting several present-day functionals of the 1-RDM. 
We compare the results with EKT, QCI(T), and experi- 
ment. 



II. THEORY 

By definition, the ionization potential and electron 
affinity are given by 

IP = E(N - 1) - E(N) 
EA = E(N) -E(N + 1), 

where E(N) is the ground-state total energy of the 
charge-neutral system, and E(N — 1), (E(N + 1)) is the 
energy of the system with one electron removed (added) . 
In the rest of the article we refer to this method of calcu- 
lating the IP and EA as the definition method (DEF). 
Due to Koopman's theorem, within the Hartree Fock 
(HF) theory, the IP in Eq. (j4|) is well approximated by 
the eigenvalue of the highest occupied molecular orbital 
(HOMO). On the other hand, within DFT, the KS en- 
ergy of the HOMO is exactly equal to the IP in Eq. Q for 
the exact xc functional. Likewise, the exact EA equals 
the orbital energy of the HOMO of the N + 1 electron 
system calculated with the exact xc functional. 
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System 




Miiller 


GU 


Power 


AC3 


PNOF1 


BBC3 


ML 


HF" 


QCI(T) 


Expt. 


He 


FIP 


25.361 


25.252 


25.579 


24.953 


24.254 


25.307 


24.626 


24.871 


24.327 


24.59'" 


H 2 


FIP 


16.490 


16.463 


16.599 


16.136 


16.419 


16.436 


16.028 


16.109 


16.245 


15.43 c 


LiH 


FIP 


8.191 


8.381 


8.490 


8.136 


8.307 


8.408 


8.027 


8.109 


7.782 


7.78 d 


H 2 


FIP 


10.259 


11.592 


11.864 


13.007 


11.614 


13.415 


12.844 


13.415 


11.919 


12.78 c 




HIP 1 


16.436 


15.619 


15.674 


16.871 


14.116 


17.361 


15.266 


15.402 




14.83 c 




HIP 2 


19.075 


17.769 


18.994 


18.721 


17.609 


19.266 


18.640 


18.857 




18.72 c 


HF 


FIP 


15.783 


15.130 


16.245 


16.925 


15.199 


17.769 


16.463 


17.116 


15.429 


16.19 c 




HIP 1 


20.381 


18.803 


20.055 


19.647 


18.591 


21.415 


20.055 


20.327 




19.90 s 


CH 4 


FIP 


13.578 


12.817 


14.123 


14.449 


12.960 


15.592 


14.395 


14.776 


14.177 


14.40 6 




HIP 1 


24.028 


23.538 


24.191 


25.497 


22.508 


20.300 


25.089 


25.633 




23.00 6 


CO2 


FIP 


9.796 


9.415 


11.320 


13.143 


10.664 


15.483 


13.633 


14.558 


13.225 


13.78 7 




HIP 1 


16.735 


16.789 


17.524 


19.102 


15.389 


19.429 


18.667 


19.075 




17.30 f 


NH 3 


FIP 


8.626 


9.878 


10.150 


10.966 


9.999 


11.510 


10.966 


11.429 


10.340 


10.80 c 




HIP 1 


17.116 


16.136 


16.626 


17.062 


15.611 


17.551 


16.599 


16.708 




16.80 6 


No 


FIP 


20.898 


20.272 


21.443 


22.585 


20.319 


23.347 


21.824 


22.640 


20.871 


21. 60 9 




HIP 1 


48.028 


47.484 


48.980 


52.600 


46.629 


52.899 


51.130 


52.219 




48.47 9 


C2H4 


FIP 


6.748 


8.218 


8.299 


9.361 


9.674 


9.796 


9.714 


10.177 


10.422 


10.68' 1 




HIP 1 


11.674 


14.068 


12.708 


14.232 


13.642 


14.340 


13.660 


13.660 




12. 80' 1 


C2H 2 


FIP 


11.048 


10.721 


10.966 


10.939 


9.821 


11.538 


11.402 


10.966 


11.184 


11.49 c 




HIP 1 


21.198 


19.701 


19.783 


19.130 


16.751 


20.136 


20.028 


18.340 




16.70 c 




HIP 2 


21.143 


22.041 


20.653 


21.633 


18.302 


20.653 


20.626 


20.789 




18.70 e 


Afip (%) 




12.38 


10.91 


7.42 


4.10 


9.19 


6.93 


2.07 


4.43 


2.91 


0.00 


Ahip (%) 




7.45 


7.29 


4.61 


8.78 


5.04 


10.83 


6.48 


6.43 




0.00 


A(%) 




10.03 


9.19 


6.08 


6.33 


7.21 


8.79 


4.17 


4.45 




0.00 






12.22 


2.47 


6.26 


2.82 


2.99 


5.16 


5.11 








A (EKT) (0/ , 




10.56 


5.69 


11.80 


6.45 


5.15 


4.21 


12.02 








A (BKT) (%) 




11.43 


4.00 


8.90 


4.55 


4.02 


4.71 


8.40 










s theorem (IP = -Bhomo) 



TABLE II: Ionization potentials (FIPs and HIPs), in eV, for various molecules calculated with different functionals using DER 
method. These results are compared with Hartree-Fock, QCI(T), and the experimental data. QCI(T) values were calculated 
with Gaussian 09 program]^ using the same basis set through Eq. Q. In the bottom row are included the percentage absolute 
average errors Afip, Ahip, and A in the calculation of the FIPs, HIPs and all IPs respectively. For comparison, the errors 
A p ™ T) , A p f IP ;T) , and A( EKT > using the EKT are also included. 



Within RDMFT, there is 110 effective single particle 
KS system reproducing the non-idempotent 1-RDM of 
the interacting system and quantities like IP and EA can 
not be obtained from an eigenvalue equation. However, 
approximate but meaningful single particle energies as- 
sociated with the natural orbitals can be defined as in 
Ref. [H]: 

e k = E({ nj })\ nh=l - E{{ nj })\ nh=0 . (5) 

The two energies on the right hand side are the energies 
of the system with all natural orbitals and occupation 
numbers having the optimal ground-state values except 



for the natural orbital of interest, k, for which occupation 
numbers are set to either = 1 or n& = 0. In this way, 
these energies are approximate electron addition and/or 
removal energies for the natural orbital k. We refer to 
this method for calculating IPs and EAs as the "energy- 
difference" method (DIF). 

It has been shown that, for extended systems [35) . the 
total energy is almost linear if a particular occupation, 
rife, is varied between zero and 1. If it was exactly linear 
then the energy difference in Eq. ([5]) would be given by 
the tangent of E({rij}). In absence of this linearity a 
good choice is to use the Slater trick and approximate e& 
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by 



<7v 



dE[{<f> 3 },{ nj }} 



dn k 



(6) 



n k =l/2 



where the derivative is calculated at the ground-state 
natural orbitals {4>j} and occupations {rij}, except for 
k which is set to rifc = 1/2. This is a good approxima- 
tion because if one expands E({rij}) at rij = 1/2 the 
term in Eq. ^ is the leading order term with the sec- 
ond order term being identically equal to zero. At first 
sight Eq. © looks similar to the Janak's theorem (43|. 
which gives the eigen-encrgies of the Kohn-Sham system 
within DFT. However, it is important to note that within 
RDMFT lack of single particle eigenvalue equations does 
not permit the direct use of Janak's theorem™ Janak's 
theorem would lead to all orbital energies, for fraction- 
ally occupied states, to be degenerate with value equal 
to the chemical potential. 



III. METHODOLOGY 



This extra variational freedom over correlates the open- 
shell systems leading to systematic errors in IPs and EAs. 
Given the exact xc functional the DEF method would be 
exact, however, for an approximate functional the DEF 
method would show these systematic errors in IPs and 
EAs. DIF and DER methods on the other hand do not 
suffer from this error since only one minimization, for 
the charge neutral system, is performed. In addition, 
it should not come as a surprise if the DER method 
performs better than DEF and DIF in many cases as 
it suffers less from possible inaccuracies introduced by 
the functional and its non unique extension to the case 
of open shells, mainly because only the 1/2 electron is 
present in the open shell. One could also consider DER 
method in conjunction with orbital relaxation (at fixed 
7i j = 1/2). However, this procedure requires full orbital 
minimization for each j making it computationally very 
demanding, while the aim of the present work is to define 
a computationally inexpensive single electron spectrum 
in terms of the optimal 1RDM of the charge neutral sys- 
tem. 



We calculate the IPs and EAs of a set of atoms and 
molecules using DEF, DIF, and DER methods [i.e., us- 
ing the Eqs. Q, (O and ©]. For comparison, results 
arc also calculated using the EKT method. Our im- 
plementation is included in a computer code for finite 
systems 44j which minimizes 1-RDM functionals with re- 
spect to occupation numbers and natural orbitals and is 
based on the expansion of the orbitals in Gaussian basis 
sets. The one- and two-electron integrals are calculated 
by use of the GAMESS program[45j|. Addition or removal 
of an electron requires the extension of the theory to o pen 
shells. In the present work, like in Refs. [l|| [30l l3ll l46| 
we use the simple extension proposed in Ref. [12J. In 
other words, we assume that orbitals arc spin indepen- 
dent while occupation numbers are spin dependent. 

For the calculation of IPs we adopt the cc-pVDZ basis 
set[47j for all elements. EAs are calculated as the IPs of 
negative ions, i.e., of N + 1 electrons. In other words, for 
both IPs and EAs, the orbital energy of the HOMO (for 
cither the neutral or ionic system) is calculated. Since 
the HOMOs of the negative ions are relatively diffuse 
states, the aug-cc-pVDZ basis set is used [ItJ . We should 
mention that a lot of neutral systems do not bind an 
extra electron. In that case, EA is equal to zero, i.e. 
the extra electron is completely delocalized. However, 
for small positive EA, the state of the extra electron can 
be delocalized and impossible to describe with localized 
basis sets. To ensure fair comparison with experiments 
a set of atoms, molecules and radicals which arc known 
experimentally to have relatively large and positive EA 
is used. 

Calculation of IPs and EAs with DEF method requires 
the difference of two energies, one for a closed-shell and 
another for an open shell-system. The broken spin sym- 
metry in open-shell systems leads to twice as many vari- 
ational parameters as there are in a closed shell system. 



Functional 


A DEF (%) 


A D i F (%) 


A DE r(%) 


Miiller 


69.85 


88.55 


63.24 


GU 


44.91 


65.96 


70.81 


Power 


64.69 


63.57 


48.73 


AC 3 


42.40 


31.17 


21.39 


PNOF1 


39.68 


61.44 


28.63 


BBC3 


30.00 


55.33 


52.00 


ML 


39.32 


47.81 


21.53 


HF 


55.00 






CI/QCI(T) 


17.28 







TABLE III: Average absolute errors Adbf, Adif, and Ader, 
in the calculation of EAs for a set of atoms, molecules and 
radicals calculated with various xc functionals using DEF, 
DIF, and DER respectively. 

Another point to be considered is that the applica- 
tion of the DER method requires the total energy func- 
tional to be continuous at rik = 1/2. However, there are 
functionals that introduce a discontinuity at = 1/2 
to distinguish between strongly and weakly occupied 
orbitals 34 1. In all cases studied here, we do not find 
optimal occupation numbers equal to 1/2. Thus, the 
step function can be safely shifted slightly away from 
rik = 1/2 without affecting the results. However, this 
procedure cannot be used in cases with optimal occupa- 
tions equal to 1/2, like H2 at the dissociation limit or 
when they vary continuously from 1 to zero. 



IV. RESULTS 

Our results for the average absolute errors in the cal- 
culation of IPs with the three methods are included in 
Table |H Average errors in the results obtained with EKT 
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FIG. 1: (Color online) The total energy, E, as a function of the occupation number, n, of the HOMO, for Ne atom and CH4 
(top) and the negative ions F~ and Cl _ (bottom). The line giving the correct experimental values for IP/EA is also shown. 
Curves are shifted to coincide at n = 1/2 in order to compare the tangents with the straight line reproducing the experimental 
results. The values at the two ends are used in DIF for the calculation of IP for Ne, CH4 and the affinity of F, CI. The 
derivatives at n = 1/2 are used for the calculation of the same quantities with the DER method. 



method are also included in the table. The actual val- 
ues for IPs obtained using the DER method are com- 
piled in Tabic HU (Full results for all methods as well as 
EKT can be found in the supplementary material [50].) 
It is clear from Table U that all functionals in combina- 
tion with the DER method give reasonable results for 
IPs with errors ranging from 4 to 13%. ML, AC3 and 
power functionals perform slightly better by giving an 
average error of only 4-6%. For the systems considered 
here, the ML functional with the DER method is the 
most accurate for the FIPs (with an error of only 2%). 
It is important to note that the errors from the DER 
and EKT method are of the same order (using the same 
functionals and basis set), while there is less computa- 
tional effort involved in DER method. The comparison 
of the DIF and DER methods allows us to assess the 
validity of the linear approximation-the linearity of the 
total energy with respect to variation of one occupation 
number while keeping the rest of the occupation numbers 
as well as natural orbitals frozen (this was demonstrated 



for solids in Rcf. [Hf). As we see in Tabic HH the DER 
method gives good results also for finite systems and for 
the best performing functionals the average difference be- 
tween the DIF and DER methods' results is in the range 
of 2-7%. This percentage may be regarded as the magni- 
tude of non-linearity of the total energy with respect to 
variation of a single occupation number. 

As mentioned in Sec. Mil the DEF method suffers 
from the over correlation error due to the approximate 
nature of the xc functional and the difference in varia- 
tional freedom between closed and open-shell systems. 
Since the DER method is less prone to this error, it is 
not a surprise that for the power and ML functionals the 
DER method improves the results over the DEF method. 
This indicates that the dependence of the total energy 
on a particular occupation number deviates from linear- 
ity but this deviation works in favor of the DER method 
by further improving the results. In order to understand 
this, in Fig. QJtop), we show the tangents at 1/2 of the 
total energy as a function of the occupation number of 
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°Ref. 4£| 
'Ref. 4C 
c Ref. 15 



System 


Miiller 


GU 


Power 


AC3 


PNOF1 


BBC3 


ML 


QCI(T) 


Expt. 


LiH 





0.192 





0.219 


0.317 


0.106 


0.264 


0.317 


0.34° 


OH 


3.157 


6.735 


3.456 


2.626 


1.784 


3.157 


2.441 


1.645 


1.83 6 


F 


5.651 


3.297 


5.040 


3.196 


3.552 


5.539 


4.363 


3.241 


3.34 6 


Li 





0.403 


0.139 


0.275 


0.400 





0.212 


0.601 


0.62 6 


CI 


3.687 


2.190 


3.632 


3.678 


2.660 


4.186 


3.623 


3.430 


3.61 6 


CN 


1.807 


2.328 


2.332 


4.264 


2.470 


4.873 


4.585 


3.627 


3.77 6 


c 2 


1.044 


1.419 


1.298 


3.386 


4.416 


4.560 


3.722 


3.055 


3.54 6 


BO 


1.243 


1.326 


1.982 


3.252 


2.076 


3.222 


3.252 


2.359 


2.83 c 


SH 


1.386 


0.706 


1.908 


1.968 


1.352 


2.454 


2.168 


2.119 


2.32'' 


PH 


0.304 





1.200 


1.252 


0.186 


2.192 


1.145 


2.023 


1.00 c 


A(%) 


63.24 


70.81 


48.73 


21.39 


28.63 


52.00 


21.53 


17.28 





TABLE IV: Electron Affinities for various atoms, molecules and radicals calculated as the IP of the system of N + 1 electrons 
with different functionals using the DER method. For systems where N + l electrons energy is found higher than the N electron 
energy, zero affinity is assumed. The results are compared with QCI(T) and experiments. QCI(T) values were calculated with 
Gaussian 09 program]!^. 



the HOMO, while keeping the rest frozen. The plots 
are made for various functionals. It is clear from Fig. Q] 
that although the total energy functional itself deviates 
from linearity the tangent at 1/2 is very close to the one 
that reproduces the experimental results. One reason for 
this improvement over the DIF method might be that 
the extension of the theory to open shells in the case the 
DER method is minimal- since only half of an electron is 
unpaired- and DER method reduces the error introduced 
by the extension of functionals to open shells. 

The average percentage errors and the values of HIPs 
for various atoms and molecules arc also presented in 
Table |T] and Table [TT| respectively. The average error in 
HIPs obtained using the Miiller, GU, power and PNOF1 
functionals, is substantially lower (4-7%) than for the 
FIPs. The rest of the functionals are less accurate for 
HIPs with average absolute errors slightly higher than 
those for the FIPs (5-10%). 

The average absolute errors in the calculation of EAs 
with the DEF, DIF, and DER methods are shown in Ta- 
ble IIII1 The actual values for EAs obtained using DER 
method are included in Table ITVl As already mentioned, 
EAs are more difficult quantities to calculate- the errors 
are introduced by describing negative ions with localized 
basis which are usually optimized for the description of 
the ground states of neutral systems. In addition, be- 
ing a small quantity, EAs are more prone to errors in 
the differences of total energies corresponding to two dif- 
ferent shell structures (for example the difference in en- 
ergy between a doublet and a singlet state). Thus, it is 
not surprising that the errors in Tables IIIII and IIVI are 
substantially higher than for the IPs. However, state- 
of-the art quantum chemical methods like QCI(T) also 
exhibit large errors within the adopted basis set. Under 
these considerations the AC3, PNOF1 and ML function- 



als perform surprisingly well for EAs with average errors 
of 21.4%, 28.6%, and 21.5% respectively. These errors 
are close to that of QCI(T) (17.3%). It is interesting to 
note that (see Tabic IIV|) there are only a few cases (5 
of 60) for which Miiller, GU, BBC3 and power function- 
als give a zero EA [E(N + 1) > E(N)], i.e., the system 
is not predicted by the corresponding functional to bind 
an extra electron. For the best performers like the AC3, 
PNOF1, and ML functionals, no such case exists. 

In order to compare the DIF and DER methods for 
the determination of EAs, Fig. [IJbottom) shows the tan- 
gents at 1/2 to the dependence of the total energy on the 
occupation of the HOMO of the negative ions F~ and 
Cl _ . The exact tangent that reproduces the experimen- 
tal EAs is also shown in the figure. Again, as in the case 
of IPs, the DER method not only looks like a reasonable 
approximation but it also improves on the results of the 
DEF method (for the functionals considered and in all 
cases studied in the present work). In particular, for the 
case of Cl~ [see Fig.QJbottom)] the tangents at 1/2 are in 
very good agreement with the exact tangent that repro- 
duces experimental EA, although the dependence of the 
total energy on the HOMO occupation number deviates 
significantly from linearity. 

A striking example of the pathological behavior men- 
tioned in Sec. IIIII is the negative ion F~. This system is 
found experimentally to be energetically lower than the 
neutral F atom by 3.34 eV (see Table I IV]) . All functionals 
underestimate this energy substantially (see supplemen- 
tary material [50j ). as a result of the enhanced variational 
freedom of the open-shell, neutral F atom compared to 
the closed-shell F~. In two extreme cases (ML and AC3 
functionals) F~ is found energetically above the neutral 
F atom. 
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V. SUMMARY 

In summary we examined the performance of the 
derivative-method proposed in Ref. [35| for calculation 
of the IPs and EAs of finite systems. The accuracy of 
IPs and EAs calculated using the derivative-method are 
compared to the IPs and EAs calculated using the defini- 
tion of these quantities [which involves two total energy 
minimizations for the system and the positive (for IP) 
or negative ion (for EA)]. In order to have a complete 
analysis we have also considered an intermediate method 
(difference method), in which the IPs and EAs are deter- 
mined by the difference of the total energies with fixing 
one occupation number to 1 and/or zero. All these results 
are further compared to the state-of-the-art CI method 
as well as experiments. 

We find that, for IPs both the difference and deriva- 
tive methods are good approximations to the definition of 



this quantity. Furthermore, it was found that the deriva- 
tive method results, obtained using Miiller, power, and 
ML functionals, are better than the values obtained using 
the definition method itself (with errors of the order of 
4-8% only). Among these functionals the ML functional 
in conjunction with derivative method is most accurate 
with errors of only up to 2%. For the EAs the errors are 
significantly larger, with the ML functional in conjunc- 
tion with the derivative method being the most accurate 
(with an error of 21%). The errors in EAs were found to 
be comparable to the errors in the CI results. 

From the present study we conclude that the results of 
the derivative method for IPs and EAs are in good agree- 
ment with experiments and this method is a promising 
technique to obtain the single-electron spectrum for sys- 
tems where state-of-the-art quantum chemical methods 
can not be applied and DFT results deviate significantly 
from experiment. 
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